Equilibria of flat and round galactic disks. 



C. Pichon 1 2 and D. Lynden-Bell 1 

1: Institute of Astronomy, The observatories, Cambridge CB30HA, UK 
2: CITA, 60 St. George Street, Toronto, Ontario M5S 1A7, Canada. 



ABSTRACT 



A general method is presented for constructing distribution functions for flat systems 
whose surface density and Toomre's Q number profile is given. The purpose of these func- 
tions is to provide plausible galactic models and assess their critical stability with respect 
to global non axi-symmetric modes. The derivation may be carried out for an azimuthal 
velocity distribution (or a given specific energy distribution) which may either be observed 
or chosen to match a specified temperature profile. Distribution functions describing sta- 
ble models with realistic velocity distributions for power law disks, the Isochrone and the 
Kuzmin disks are provided. Specially simple inversion formulae are also given for finding 
distribution functions for flat systems whose surface densities are known. 
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1. Introduction 

Over the next decade, ample and accurate observational data on the detailed kinematics of nearby disk 
galaxies will become available. It will be of great interest to link these observations with theoretical models 
for the underlying dynamics. The problem of finding the distribution function for an axially symmetrical 
system may formally be solved by using Laplace transforms (Lynden-Bell 1960 [jT2| ) or by using power 
series (Fricke 1952 ||, or more recently Vauterin 1995 ||) . Kalnajs (1976) Jll[, followed by Miyamoto 
(1974) jl4| chose specific forms of distribution function because the flat problem has no unique solution. 
Distribution functions for the power-law disks were also recently derived by Evans (1993) while Hunter & 
Qian (1993) || presented an inversion scheme which can be applied for thickened disks with two integrals. 
For flat disks, it is desirable to use the functional freedom left in / in order to construct a distribution 
function which accounts for all the kinematics, either observed or desired (i.e. which accounts for the line 
profiles observed, or which are marginally stable to radial modes). Independent measurement of the observed 
radial and azimuthal velocity distribution functions could, for instance, be contrasted with predictions arising 
from the gravitational nature of the interaction. Indeed the laws of the motion and the associated conserved 
quantities together with the assumption that the system is stationary put strong constraints on the possible 
velocity distributions. This is formally expressed by the existence of an underlying distribution function 
which characterises the dynamics completely. The determination of realistic distribution functions which 
could account for observed line profiles is therefore an important project vis a vis the understanding of 
galactic structure. Producing theoretical models that accounts globally for the observed line profiles of a 
given galaxy provides a unique opportunity to inspect the current understanding of the dynamics of SO 
galaxies. It should then be possible to study quantitatively all departures from the flat axisymmetric stellar 
models. Indeed, axisymmetric distribution functions are the building blocks of all sophisticated stability 
analyses, and a good phase space portrait of the unperturbed configuration is often needed in order to asses 
the stability of a given equilibrium state. Numerical N-body simulations also require sets of initial conditions 
which should reflect the nature of the equilibrium. 

For a flat galaxy all the stellar orbits are confined to a plane and by Jeans' theorem the steady state 
mass- weighted distribution function must be of the form / = /(e, h), where the specific energy, e, and the 
specific angular momentum, h, are given by 




(1.1) 
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The surface density, arises from this distribution of stellar orbits provided that the integral of / overall 

bound velocities is i.e. 

S(iJ)= / / f{e,h)dv R dv^, (1.2) 



where the integral is over the region |(vfj < V- The following inversion methods assume that S(i?) is 

given and its potential on the plane, ip(R) is known. This is achieved either by requiring self-consistency via 
Poisson's equation 

V 2 ^ = -4tt(5(z) E(R) , (1.3) 
or alternatively by assuming that the disk is embedded in a halo and choosing both £ and ip independently. 



At constant v§ and R, v R dv R — de and hence, using Eq. (1.1) 



dv R = de[2(e + Tp)-h 2 R-' 2 ] 1/2 . (1.4) 



Furthermore, at constant R, dh = Rdv^, and therefore [Eq. (1.2)| yields 



f +^/2R^p ,o ft h) dedh 

Z(R)= 2/ J[e,n) aean — 

The factor of 2 arises because Vi? takes both positive and negative values which are mapped into the same 
range of e which depends only on v \ . Similarly, 



2 / (e, h) y/2 (e + -0) i? 2 — /i 2 de d/i . (1.6) 

The b asic inver sion problem is concerned with finding the class of f(e, h) compatible with the constraints 



1.5)| and |(1.6) 



Two classes of inversion methods are constructed and implemented in this paper. These lead to classical 
distribution functions compatible with a given a given surface density (section 3), or a given surface density 
and a given pressure profile (section 2). The former technique relies on an Ansatz and yield direct and 
general methods for the construction of distribution functions for the purpose of theoretical modelling. The 
latter technique yields distribution functions that correspond to given line profiles (or rather, to the shape 
of the line profiles induced by the velocity distributions which will here loosely be called line profiles) which 
may either be postulated or chosen to match the observations. In fact, only a sub-sample of the observation 
is required to construct a fully self-consistent model which accounts, in turn, for all the observed kinematics. 
The corresponding redundancy in the data may be exploited (as sketched in section 4) to address the 
limitations of the description. 



2. Distribution functions for given kinematics. 

This global inversion method introduces an intermediate observable, F^ > (R,v ( p), the distribution function for 
the number of stars which have azimuthal velocity at radius R (or alternatively G £ (e,R) the distribution 
function for the stars which have specific energy e at radius R). This technique may be applied to the 
reconstruction of distribution functions accounting for observed line profiles. In fact, given some symmetry 
assumptions about the shape of an observed galaxy, it is shown that only a subset of the available line profiles 
- say the azimuthal velocity distribution - is required to re-derive the complete kinematics. Therefore this 
method provides a general procedure for constructing self-consistent models for the complete dynamics of 
disk galaxies. These models predict radial velocity distributions which may, in turn, be compared with 
observations as discussed in section 4. 

Alternatively, a 'natural' functional form for or G £ may be postulated and parameterised so that it is 
compatible with imposing the surface density, the average azimuthal velocity and the azimuthal pressure 
profiles (or equivalently the Toomre number Q, as Q follows from the equation of radial support). The 
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distribution function f(e, h) of this disk follows in turn from or G e via simple Abel transforms. This 
prescription is very general and especially useful when setting the initial conditions of numerical N-body 
stability analysis. 



2.1. Inversion via line profiles 

The number of stars which have azimuthal velocity within dvj, at radius R reads 



F <t ,{R,v <t> ) = J f(e,h)dv R =V2 J 



f(e,h) 



de , 



(2.1) 



where the effective potential, Y, is given by Y — ip — h 2 R 2 /2 = ip — w|/2. The line profile F ( f,(R,v l j > ) 
may be expressed in terms of (h,Y). Indeed the identity Y = if) — h 2 R~ 2 /2 may be solved for R which 
yields R(h,Y). The azimuthal velocity v$ becomes in turn a function of h, Y defined b y = h /R(h,Y) 
for each branch corresponding to the two roots of Y = ip — h 2 R~ 2 /2 (as illustrated in Fig. 2.1). Calling 
F$(h, Y) = F C / ) (R, v,/,) allows the inversion of Eq. (2.1) by an Abel transform: 



1 



dF* 



dY 



V2ir J dY (_ £ ) _ Y 



(2.2) 



where the partial derivative is taken at constant h. It was assumed here that the distribution F^R^v^) 
vanishes at the escape velocity. Note that F yields both the symmetric and the antisymmetric parts of the 
distribution function. The r.h.s of Eq. (2.2) is advantageously re-expressed in terms of the variables (R, h) 
given that (for monotonic integration) 



dF^ 
dY 



dY 



dFj, 
OR 



dR 
dY 



dY 



dFj, 
dR 



sign 



dR 
&Y 



dR. 



(2.3) 



This yields two contributions for Eq. (2.2 

R 



{dF <t> /dR) h dR 
^/h 2 /R 2 -2ij(R)~2e 



(dF <t> /dR) h dR 
y/h 2 /R 2 -2ij(R)-2e 



(2.4) 



where R p (h,e), and R a {h,e) are, respectively, the apogee and perigee of the star with invariants h and e, 
and R e (h) is the inner radius of a star on a "parabolic" (zero energy) orbit with momentum h. Note that 
the derivative in Eq. (2.4) is performed keeping h = Rv^ constant. The contributions to the two integrals is 



illustrated on Fig. 2.1 



Equation (2.4) yields the unique distribution function compatible with an observed azimuthal velocity dis- 
tribution. All macroscopic properties of the flow follow from /. For instance the surface density reads 



E = 



/+ (e,h)dv R 



dvrfj — 2 J Frf, (R, vj,) dv<f, , 



(2.5) 



while the radial velocity distribution, Fr(R,vr), obeys 



Fr (R, vr) 



1 

R 



f 



h 2 
2R 2 



ip(R),h) dh 



(2.6) 



-R(2i,-v R ) 



Both Fr and E may in turn be compared with the corresponding observed quantities. 
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Figure 2.1: the relationship between integration in terms of the effective potential, Y, 
and R integration. The effective potential is drawn here as a function of R for two values 
of h. The area between that curve and the lines Y = — e and R = lZ e is shaded, defining 
3 regions from left to right. The middle region does not contribute to Eq. (2.4). The 
equation Y = — e has two roots corresponding to the perigee and the apogee of the star, 
while Y — has two roots corresponding to infinity and lZ e , the inner bound of an orbit 
with zero energy and angular momentum h. The sign of the slope of Y(R) gives the sign 
of the contribution for each branch in Eq. (2.4). 

2.2. Disks with given Q profiles 

From the point of view of stability analysis, F$ may be chosen to match given constraints such as a specified 
potential and temperature profile of a disk. The azimuthal pressure, p^ = — (i^) 2 ), is then fixed via 

the equation of radial support and by the temperature of the disk defined by Toomre's Q number: 



d(R PR ) 
OR 



2v3 



and pr = 



7T 2 £ 



where the circular velocity, the epicyclic frequency and the surface density are given by 

\2" 



v 2 = 



1 dtp 
RdR 



1 
if 1 



d{RV c y 



OR 



and S = 



1 

27 



dtfj 
dz 



(2.7) 



(2.8) 



Note that here the azimuthal stress p^, contains the mean streaming stress. It is assumed here that the 
field is self-consistent and obeys Eq. (1.3); alternatively S may be given independently of tp and the method 



described here still applies. The surface density, S, and the azimuthal pressure, £( 
expressed in terms of F ( j,(R, v^) via 



may in turn be 



/+ (e,h)dv R 



dvA 



F</> (R, dvj, , 



(2.9a) 
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/+ {e,h)dv R 



,dv<p 



F<j, (R, v$) v\ dvcf, ■ 



(2.96) 



where the circular escape velocity, v e , is equal to y/2yj. The function / + stands for the component of the 
distribution function even in h. Any function satisfying these moment equations corresponds to a state 
of equilibrium stable against ring formation when Q > 1. Realistic choices for are presented in the next 
sections. 



2.3. Alternative inversion method 



Another intermediate observable, G e (i?,e), the distribution function for the number of stars which have 
specific energy e, may also be parameterised to fix the surface density and the average energy density 
profiles. For external galaxies G e is not directly observed. For our Galaxy, G £ is measurable indirectly - via 
the distribution of stars with given radial and azimuthal velocity given by spectroscopy and proper motions 
- but only yields the even component in h of the distribution function. However, the corresponding inversion 
method still provides a route for the construction of models with specified temperature profiles. 
The number of stars which have specific energy s within de at radius R reads 



G e {R,e)= f!±& 

J V R 



f+(e,h) 2 
dv v = — 



g(g,fe 2 /2) 
y/2 J y/X - h 2 /2 



where X is defined by 



X = R 2 (if) + e) . 
Here the auxiliary function g is given by 

g(e,h 2 /2)=U(e,h)/\h\. 

The surface density, and mean energy density, (e), may be expressed in terms of G e (R,e) via 

o 



S = 



/+ (e,h)dv tj) dv R = J G e {R,e)de, 



S (e) 



f+{e,h) 



2 



dv<j> dv R 



G e (R, e) ede . 



(2.10) 

(2.11) 
(2.12) 

(2.13a) 
(2.136) 



Note that only the symmetric component of the distribution function follows from Eq. (2.12). The local 
mean energy density (e) = S _1 (p/j/2 + — ip{R) is fixed by the temperature of the disk defined by 

Toomre's Q number: 



d(Rp R ) 
OR 



PR 



- ip (R) given p R , = Q : 



2v3 



7T S 



(2.14) 



where n and V c are given by Eq. (2.8). The function G(R, e) may be expressed in terms of e, X via Eq. (2.11) 
and R — R(ip). Calling G e (R,X) = G e (R,e) leads to the inversion of Eq. (2.10) by an Abel transform: 



h 2 /2 



g(e,h 2 /2) 



dG e 



dX 



V2ir J dX y /h 2 /2-X 



(2.15) 



/+(e, h) follows from g(s, h 2 /2) via Eq. (2.12). 

Any function G e satisfying the moment constraint Eqs. (2.13) yields, via Eq. (2.15) and (2.14), a maximally 
rotating disk stable against ring formation. 
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2.4. Implementation: Gaussian line profiles 

Gaussian velocity distributions are desirable both as building blocks to fit measured line profiles and as 
"realistic" choices for F$ in the construction scheme of disks parameterised by their temperature. The 
former point is discussed in section 4. The construction of Gaussian line profiles compatible with a given 
temperature requires a supplementary assumption for the mean azimuthal velocity of the flow, (tu), on 
which the Gaussian should be centred. In additions to the the two constraints, Eqs. (2.9), this puts a third 
constraint on i 7 ^, namely 



/_ {e,h)dv t 



v<$, dv^ — 2 J F$ (R, ity) v<t, dv<p . 




(2.16) 



where /_ is the odd component in h of the distribution function. Suppose the following functional form for 

l2\ 



F^ (R, tj ) = S (R) W n (R, V+) exp 



[v$ - v (R) } 
2a 2 (R) 



(2.17) 



where the window function W n is intended to damp F c / ) (R,v ( p) near the escape velocity v e . A possible 
expression for W n is 



W n {R,V^) = 



exp 



(yl (R) v; 



if M < V e , 



elsewhere . 



(2.18) 



Self-consistency requires that the unknown functions (v, a, S) of Eq. (2.17) are solved in terms of (v^) , S, and ; 
via Eqs. (2.9) and Eq. (2.16")} For practical purposes, the temperature range explored in realistic disk models 
is such that (v e — v) 2 /2a 2 is quite large at all radii; if n is also chosen so that at temperature a 



n > 



(2.19) 



then the tail of the Gaussian function need not be taken explicitly into account. The line profile F then 
reads directly in terms of (v^) , S, and p^: 



(R, v 4 ) 



\/27r <T0 



exp 



2a 2 



(2.20) 



where o§ is the the azimuthal velocity dispersion 



(2.21) 



The azimuthal pressure follows from the equation of radial support and the temperature of the disk 
defined by Toomre's Q number given by Eq. (2.7) . The expression of the average azimuthal velocity, (u^), 
may be taken to be that which leads to an exact asymmetric drift equation: 



K 2 R 2 



(2.22) 



Equations (2.20)-(2.22), together with Eqs. (2.7)1 provide a prescription for the Gaussian azimuthal line 
profile which yield the distribution function f(e, h) of a disk at given temperature profile via Eq. (2.4). 
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Figure 2.2: isocontours of the Gaussian distribution functions for the Kuzmin disk with 
temperature: Q = 1, (top left panel), Q — 1.25, (top right panel), Q = 1.75 (bottom left 
panel ), Q = 2 (bottom right panel). The construction scheme is described by Eq. (2.4) 
and Eqs. (2.23). Each diagram represents the number of star with angular momentum, 
h, and relative energy, r/ = e/eh, where Sh is the energy of the star on a circular orbit 
with momentum h. The parametre rj therefore measures the eccentricity of the orbits. A 
hotter component is apparent at larger momentum for the Q = 2 disk where the contours 
are more widely spaced. 



2.5. Example: constant temperature Gaussian Kuzmin disks 

For the Kuzmin disk, the prescription described in the above for the parameters of the line profile yields 



O"0 = 



Q 



4(1 + i? 2 ) 



3/4 



and (v<p) 



R (256 - 60 Q 2 + 128 R 2 - 21 Q 2 R 2 + 16 R 4 ) 
4(l + i? 2 ) 3/4 (4 + i? 2 ) 



1/2 
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where Q is Toomre's number. From Eq. (2.17), F$ becomes as a function of h 



[h, R] = 



VW 3 



(4 + R 2 ) 



-exp 



8(l + i? 2 ) 3/2 ( h R (256 - 60 Q 2 + 128 R 2 - 21 Q 2 R 2 + 16 R 4 ) 1/2 ' 



Q 2 



Differentiating with respect to R gives 



d log Fa 



dR 



(4 + R 2 ) 



R 



8 h 2 (2 - i? 2 ) 
Q 2 i? 3 (l + i? 2 ) 1/4 



4(l + i? 2 ) 3/4 (4 + i? 2 ) 



ARP t 



? 2 (l + i? 2 ) 3/4 (4 + i? 2 ) 3 
6hF 6 R 



with 



2(l+i? 2 ) 7/4 Q 2 (1 + i? 2 ) (4 + i? 2 ) 
(-1024 + 216 Q 2 - 768 i? 2 + 106 Q 2 R 2 - 192 R i + 7Q 2 R i - 16 i? 6 ) 



(256 - 60 Q 2 + 128 R 2 ~ 21 Q 2 i? 2 + 16 R 4 ) 



1/2 



and P 6 = -256 + 60 - 192 R z + 27Q 2 i? 2 - 48 i? 4 - 4 i? 1 



(2.23) 



(2.24a) 
(2.246) 



The equations (2.23) and (2.24) together with Eq. (2.9) fully characteri se a com plete family of distribu- 
tion functions parameterised by their temperature via Toomre's num ber. Fig. 2.2 illustrates this inversion 
procedure. The parametrized azimuthal velocity profile is given on Fig. 2.3, together with their derived 
radial velocity profile. In contrast, Fig. 2.4 corresponds to the same method applied to the Isochrone disk, 
ijj = 1/(1 + y(l + R 2 ))- Note that the relatively broader core of the Isochrone disk does not require as strong 
a hot component to match the imposed temperature profile. 




Figure 2.3: Azimuthal (left) and radial (right) velocity distributions for the Gaussian 
Kuzmin-Toomre disks parametrized by their "Q" temperature. The left panel was pa- 
rameterised and the right deduced using Eq. (2.6) and Eq. (2.4). The curves from top to 
bottom correspond to radii going from R = 0.5, 1, • • • 2.0. The plain curve corresponds to 
a Q = 1 disk, the dashed curve to a Q = 1.25 disk, while the doted curves corresponds 
to a Q = 1.5 disk. Note that, as expected, the hotter disks have broader radial velocity 
distributions. 
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Figure 2.4: isocontours of the Gaussian distribution functions for the Isochrone disk 



displayed as in Fig. 2.2 



2.6. Example: Power law disks 

The potential and surface density for the power law disk read 



iP = R~? and V=2£ R -P-i, 
2tt 



(0 < p < 1) 



where Sp is given by : 



/ r[l/2 + /3/2] T[l-{3/2} \ 



(2.25) 



(2.26) 



It is assumed here that distances are expressed in terms of Rq, the reference radius, and energies in terms of 
ipo — i/j(Rq). (in units of G = 1). The pressures follow from Eq. (2.25) and the equation of radial support 
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Eg. (2.7)} The pressures follow from Eq. (2.25) and the equation of radial support 

Q 2 sl 



PR 



8tt/3 (2 - /3) 



R 



-2/3-1 



2tt 



Q 2 s 2 

2/3(2-/3) 



R 



2/3-1 



where Q = (t r k/(ttGY,). Assuming that F^h, Y) scales like CY n+1 h b , the inversion yields 



f(e,h) 



dF 6 



dY 



(n + 1) CT (n + 1) ( _ £) n+i/2 hb 



V2ir J dY y/(-e)-Y v^) r (n + 3/2) 
The surface density, E, and the azimuthal pressure, S(w?), read in terms of F^Qi, Y) via 

E = 2 /' A,, = g2^r(2 + n)r(6/2 + l/2) fi _, (3/2+ „ +fe/2) - 6 
* * (n + 1) v^r)r (6/2 + n + 5/2) 



(2.27) 



(2.28) 



o 



E(^\ =2 Fsvldv 



C2 fe/2 + 5/2 r(2 + n)r(b/2 + 3/2) fl _ /3(5/2+w+fe/2) _ 6 _ 1 



(n + 1) V(2^) r ( fo /2 + n + 7/2) 



(2.29a) 



(2.296) 



where the circular escape velocity, v e , is equal to \/2nJ>. Identifying Eqs. (2.29) with Eqs. (2.25)| -(2.27) yields 
the relations 



b = 



2(3n 



1, n = /3-2 + 



2(/3-2) 2 /3 



2-/3 

which gives the distribution function 



Q 2 sl 



c = 



Sp2-l 3 ' 2 - 1 l 2 T(n + /3/2 + 5/2) 



f + {e,h)=Sp (- £ ) n+1/2 |A| 



27rr (2 + n) T (/3/2 + 1/2) 

2 -5/2 ^-3/2 p [Q , + 



r[n + 3/2]T[a-n - 1] 



2\ /3n/(2-/3)-l 



(2.30) 



(2.31) 



where a = 1 4- 2n/(2 — /3). Here, the effect of the parameter n on temperature becomes clear when inverting 
Eq. (2.30) for Q: 

1/2 

(2.32) 



2-/3 



2/3 



_2-/3 + ?i_ 

These results for the power-law disks were also derived by Evans H who used the Ansatz II described in 
section 3.1 and illustrated in section 3.3. 



3. Distribution functions for a given surface density 

As mentioned in section 1, the inversion problem for flattened disks has no unique solutions. Here the 
freedom in choosing the distribution function is exploited to choose special forms of distribution function 
which make the solution of the resulting integral equation remarkably simple. Explicit inversions are given 
for distribution functions whose even parts are of the form f + (e,h) — (— e) n+ 2 G n (h 2 /2) where h is the 
specific angular momentum and e the specific energy of a given star. Inversions are also provided when 
f+(e, h) = h 2n F n (e). Examples are given and illustrated graphically . The part of / which is antisymmetric 
in h = Rv,), cannot be determined from the surface density alone unless it is assumed that all the stars rotate 
in the same sense about the galaxy, in which case /_ = sign(ft,)/ + . More generally, if the mean velocity 
(Vtj,) is specified - say by the asymmetric drift prescription - a similar Ansatz for /_ allows its explicit 
determination too. The inversion formulae require expressing powers of R times the surface density E(i?) 
as functions of either the potential ip(R) or of the related function Z(R) = R 2 ip. These formulae can only 
be applied to those density distributions for which the potential in the plane of the matter is known (either 
self-consistently or not). The other characteristics of the disk follow from the knowledge of / and cannot be 
specified a-priori. 
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3.1. Derivation 



• Even component of f: first Ansatz 

A first Ansatz is to look for solutions for / + which take the form 



f + (e,h) = (-s) n+ i\h\G n 



(3.1) 



where n is first assumed to be an integer. Here the \h\ is taken out for convenience; its appearance in no way 
implies that /+(e, 0) = because G n (h 2 /2) may behave as for small h. The quantity n measures the 

level of anisotropy in the disk. Large n correspond to cold, centrifugally supported disks; small n correspond 
to isotropic, pressure supported disks. Inserting the Ansatz into Eq. (1.5) gives 



2i 
R 



R 2 i> 



G r , 



y/Y-(-e) 



h 2 



(3.2) 



The inner integral can be transformed to Y n+l I n where 



o 7t 



2 dx 



^r(n + |) _ (2n+ 1) (2n- 1) ...1 



2r (n + 3) 



2 n+2 ( n+ i)i 



(3.3) 



with x — (— e/Y). Re-expressing Y n+1 in terms of the potential, writing H for h 2 /2 and Z for R 2 ip, this 
expression becomes on multiplication by R 2n + 3 



S n = i? 2n+3 S (R) = 2§/„ / (Z - H) n+1 G n (H) dH , 

Jo 



(3.4) 



hence defi ning S n . Z = R 2 ij} is known, so R 2n +3 I](R) can be re-expressed as a function S n (Z). Differentiating 
n + 2 times with respect to Z , using Eq. (3.3) and substituting H for Z gives: 



Eq. (3.4) 



G n (H) 



1 



V2[(r 



d 

dk 



^ Sn(H) 



(3.5) 



Eq. (3.5)| specifies the part of the distribution function which is even in h. Demanding no counter-rotating 
stars would imply, for instance, / = / + + /_ =0 for h < 0. In that case /_ = — / + for h < which leads to 
/_ = sign(h) f + (e,h). These solutions are called maximally rotating disks. The formal solution, Eq. (3.5)| , 
must be non-negative for all h so that Eq. (1.5) corresponds to a realizable distribution of stars. Within that 



constraint, the choice of n is, in principle, free. There are many distribution functions which give the same 
surface density distribution because there are two functional freedoms in f(e, h) and only one function S(i?) 
is given. Here it has been shown that even within the special functional form of Ansat z I there i s potentially 
a solution for any chosen integer n. When n is not an integer, n — no — a, < a < 1, Eq. (3.2) may still be 
inverted using Abel transforms, though the final solution involves an integral which might require numerical 
evaluation 



G n (H) 



sin (na) 2- 3 / 2 tt- 1 /- 1 
[(n + 1 - a) (n - a) ... (1 - a)} 



H 



«o+2 



dZ 



dZ Sn J(H-Zf-" 



+ 



no + l 



Sr. 



1 



dZ ~V z=0 Hi-° 



(3.6) 



where /„ is defined for non integer n by the first identity in Eq. (3.12 



Alternatively, there are many ways in 
which it is possible to re-express S(i?) in the general form S(i?) = R~ 2n ~ 3 S*(R 2 4>) ■ For any particular 
sum, a distribution function G* n (H) can be found which gives rise to that part of the density given by the 



n term S*. Indeed the relationship between G* and S* is just that given by Eq. (3.5) between G n and S n 



Since the original integral equation is linear, it follows that each expression for S(i?) yields a distribution 
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function in the form f + (e,h) — J2 n (~ £ ) n+ ^ l^l^n (^V^)- There is no requirement that each G* should be 
positive provided that the sum / + is positive for all e and h. Thus these expression for the density and 



distribution functions gives a very considerable extension of the original Ansatz Eq. (3.1) 
• Odd Component of f and asymmetric drift 

When the average azimuthal velocity field is known, a similar inversion procedure is available in order 
to specify the odd part in h of the distribution function. This field may be assumed to be given by the 
asymmetric drift equation which takes the form 



£ (V = p<t> - pr [ 



(3.7) 



where the azimuthal pressure, p$, and the radial pressure, pr, are derived from the even component of /_, 
and the epicyclic frequency and the circular velocity k and V c follow from ip. Formally, the only difference 
from the previous analysis is that in place of Eq. (1.2), the integral equation relating /_ and (i^) becomes 



BE = 4 



Replacing Ansatz I by Ansatz I': 



I' 



f- (e, h) 



^2 (e + ^)R 2 -h 2 



zhde dh . 



(3.8) 



/_ (e, h) = (-e) 5 G « 1 Y ) si § n W > 



(3.9) 



yields a solution for G n which is formally equivalent to Eq. (3.5) , but with S n defined by S n (Z) = 
R 2n+i I] (v<f,). The extension to solutions for /_ when £ (v^) is a linear combination of i?^ 2n ~ 4 5*(Z) is 
straightforward. 

• Alternative Ansatz 

In place of Ansatz I, consider a different Ansatz: 



II 



f + (e,h) = h 2n F n (s) , 



(3.10) 



where n is also first assumed to be an integer. Inserting Eq. (3.10) into Eq. (1.5), and reversing the 
order of the integrations, yields 



S(J?)=4 f F n (e) f 

J —ih JO 



rl/2 



h 2n dh 



-iji .hi y/X — h 2 

where X was defined in Eq. (2.11). The inner integral is X n J„, where 



de , 



(3.11) 



1 x 2n dx v^ r (i 
o \/l - X 2 



2T(l + n) 



= [(»-s) (»-D-3] *l ( 2n! ) > 



(3.12) 



with x = h/X*. Thus [Eq. (3.11)| becomes, using Eq. (2.11) for X n 

£ (i?) = 2 n+2 J n R 2n f° (e + j,) n F n (e) . 



(3.13) 



Re-expressing R S(_R) = S n (ip) and writing e for — ip yields 



Fn(s) 



[(n-|) (n -§)..!]* Vd(-e) 



n+l 



(3.14) 
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Equation [Eg. (3.14) was derived independ ently by Sa wamura (1987) JlTy . 

When n is not an integer, the solution to Eq. (3.13) can still be found, though with a supplementary 
integration (with n = n Q — a, n the closest upper integer), to give 



Fn (s) 



-n-2 



n J„ sin net 



[(n + 1 - a) (n - a) ... (1 - a)] 



d n ° +1 

dip 
1 f d n ° 



dip 



dip 



■S, 



[e-ip] 



(3.15) 



where J n is defined for non integer n by the first identity in Eq. (3.12). The Ansatz in Eq. (3.10) may 
be also generalised when E(i?) is expressed as linear combinations of i? 2 ™6> *(?/>). 
• Alternative Ansatz: odd component 

The above procedure is also straightforward to implement in order to constrain the antisymmetric 
component of the distribution function. Indeed Ansatz II': 



II' f-(e,h) = h 2n - 1 F n (e) 



(3.16) 



may be chosen in place of Ansatz I'. Then Eq. (3.11) becomes 



EE(^)=4/ F n (e) 



h 2n dh 
VX-h 2 ' 



:de . 



(3.17) 



Hence the solution for F n is formally identical to Eq. (3.14)] but with S n substituted for <5>„, where 



Sn(ip) = R 1 2 ™S (v^). Again, the integral equation is linear, so the prescription defined by Eq. (3.14) 
when £ (v^) can expressed as linear combinations of R 2n ~ 1 S*(ip). 



3.2. Disk properties 



• Radial velocity dispersion 



The average radial velocity dispersion, <r 2 R , defined by Eq. (1.6) is an important quantity for the local 
stability of the disk. Given Eq. (3.1)| and using Ansatz I, Eq. (1.6) may be rearranged as 



2 2 2 T n 

L-iO I " 



R 2n+5 



(Z - H) n+2 G n (H) dH . 



(3.18) 



The calculation o f cr\{R) is the refore str aightforward once the function G n has been found. Note the 
similarity between Eq. (3.18) and Eq. (3.4) . In fact this similarity provides a check for self-consistency since 
one should have 

^ {a 2 R 1ZR 2n+b ) =2(n + 2)2 n //, i i? 2n+3 £ = i? 2n+3 S . (3.19) 
Putting |Eq. (3.10)| into [Eq. (1.6J yields, for Ansatz II: 



Sa| = 2 n+3 J n R 2n 



(e + ip) n+1 F n (e) de . 



(3.20) 



The similarity between Eqs. (3.13) and (3.20) yields the identity: 

d 



dip 



(Ea 2 i?- 2 ") 



2 (n + 1) Y,R- In J n /J n = Y.R 



-2 n 



(3.21) 



Azimuthal velocity dispersion 
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The average azimuthal velocity dispersion, er^, is an observable constraint on models which satisfies 

(3.22) 



(4 + M 2 ) = e (4) = ^ 



/+ (e, h) h 2 de dh 



2J?2 



2 (e + ^) - £ 



where the mean azimuthal velocity (tty) is given in terms of /_ by linear combinations of Eqs. (3.8) and 
(3.17)| . Putting |Eq. (3.1)| into Eq. (3.22) and following the substitutions |Eq. (3.1)||3~5)| yields, for the first 



Ansatz, 



R 2n+5 S ^2 + {vrf>) 2^ = 2 | ^ J* (z _ H) n+l dH 



Putting Eq. (3.10) into Eq. (3.22) gives, for the second Ansatz: 



E (v 2 ) = 2»+ 4 J n+1J R 2 " f 

J-th 



(e + ^) n+1 F n (e) de . 



(3.23) 



(3.24) 



Note that <tr, -tp, £ and (t> 2 ) are related via the equation of radial support Eq. (2.7% The properties of 
the disks constructed with Ansatz I & II are illustrated in Figs. 3.2 and Fig. 3.1 for the examples described 
below. 



3.3. Examples 



• Isochrone disks 

Distribution functions for the flat Isochrone are simplest with the assumptions of Ansatz II. In the 
equatorial plane, the Isochrone potential reads 

tp = GMj (b + r b ) where r b 2 = R 2 + b 2 . (3.25) 

The corresponding surface density is 



Mb 
2vri? 3 



{In [(R + n ) /b] - R/n} 



Mb 
2nR 3 ' 



(3.26) 



But R 2 ip = Z = (rb 2 — b 2 )GM / (r b + b) = GMb(s — 1), where s is the dimensionless variable r b /b. Therefore 



T,R 



2n+3 = ^R 2n L (s) = (s 2 - 1)" L (s) , 



(3.27) 



where L(s) — ln(Vs 2 — 1 + s) + \/s 2 — 1/s. Note that L'(s) is quite simple: L'(s) = Vs 2 — 1/s 2 . With this 
notation, 



d 



n+2 



which implies, in turn, that 



G,, 



J2 ) ( i?2 " +3E ) 



Mb 2n+1 / Q 

2tt (Gmb) n+2 \~d~s 



n+2 



(s 2 -l) n L(s) 



M 



-n—1 Ln — 1 



d 



2 5 / 2 ttG7„ (n + 1) \ds 



n+2 



[( S 2 ~l)L(s)] 



s=l+h 2 /2 

Therefore, the final distribution function (or rather its symmetric part, /+) reads 

<-e/(GM/b)} n+1/2 ' ^ 2 x 1/2 



f+(e,h) 



h 2 y^fd^ n+2 



(s 2 -l) n L(s) 



(3.28) 



(3.29) 



(3.30) 



4TT 2 Gb(n+ 1/2)!! \2GMb J \ds 

with s = 1 + h 2 /(2GMb). Using these distributions functions, Toomre's Q criterion for radial stability 
(Q = ct h k/3.36GEo) may be calculated via Eq. (3.18). Note that for the Isochrone potential 

n+2 



4 (R) 



2 5 / 2 l 



R2n+5<Z (#) J Q 



(Z - H 



n+2 



(I) 



dH. 



(3.31) 



s=H+l 




• Toomre-Kuzmin disks 

Ansatz II is more appropriate when seeking distribution functions for the Toomre-Kuzmin disks. In 
the equatorial plane, the potential of the disk reads tp = —GM/rb, while the surface density is £ = 
(2n)^ 1 GM bjr\. Ansatz II yields the distribution (in units of b and GM/b ) 



f+(e,h) 



2" 



-2 u2n 



„2n+3 



(1 



(n+1) 



(3.32) 



This example illustrates the drawbacks of Ansatz II. The factor h in [Eq. (3.32)| removes all zero angular 



momentum orbits. Self-consistency yields a solution with circular orbits near the angular momentum origin 
(/ scales like powers of h 2 /[l — (— e) 2 ] 3//2 ). In position space, this implies that the inner core of the galaxy is 
made of circular orbits! Realistic distribution functions will therefore require superpositions of solutions of 
the family defined by (3.32) , as illustrated by Fig. 3.2 . It can be shown from Eqs. (3.32) and Eq. (3.24)| that 
£ u 2 R does not fall off fast enough to yield cold disks (it ought to scale like tp 6 to lead to constant Q profiles). 
• Power law disks 

The simplicity of power law disks leads to solutions for the inversion problem which may be extended 
to models with a continuous free "temperature" parameter as described above. This temperature may be 
varied continuously between the two extremes corresponding to the isotropic and cold disk. The potential and 
surface density for the power law disk where given by Eqs. (2.25)-(2.26). The inversion method corresponding 
to Ansatz I may be carried out and, using (3.6), yields the distribution function found in Eq. (2.31). The 
velocity dispersions are recovered from Eqs. (3.18)| and (3.22): 



(3.33) 
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0.025 - 



0.02 



0.015 - 



0.01 - 



. 005 




Figure 3.2: radial (bottom lines) and azimuthal (top lines) pressure for a Kuzmin 

The weights are (2/3, 1/12, 1/4- 



disk constructed by linear superposition of Eqs. (3.32) 



i/12^/24, i/24) where i = 10, 13 • • • 22 for the n = 0, 1 
superposition has positive mass everywhere. 



models. It was checked that the 



By construction, these dispersions satisfy the equation of radial equilibrium: 



(3 (2 - (3) 
(2-p + n) 



R- 



d (RT,a 2 R ) 



(3.34) 



4. Observational implications 

The inversion method described in section 2 may be applied on observational data which are now becoming 
available (e.g. Bender et al. (1994) [jlj and Fisher (1994) Q have presented line of sight velocity distribu- 
tions of elliptical and SO galaxies obtained from spectra with (60 km/s) resolution and S/N ~ 50 ). The 
observations may be carried out as follows: consider a disk galaxy seen almost edge on. It should be chosen 
so that it looks approximately axisymmetric and flat. It may contain a fraction of gas in order to derive the 
mean potential ip from the observed H I velocity curve (Alternatively, the potential may be derived directly 
from the kinematics if the asymmetric drift assumption holds). Putting a slit on its major axis and cross 
correlating the derived absorption lines with template stellar spectra yields an estimate of the projected ve- 
locity distribution which should essentially correspond to the azimuthal line profile F^(R,va). An estimate 
of the radial line profile Fr(R,Vr) arises when the slit is placed along the minor axis. But Fr(R,vr), is 
given by Eq. (2.6) while Eq. (2.4) gave / expressed as a function of F^v^, R). Therefore, putting Eq. (2.4) 
into Eq. (2.6) provides a simple way of predicting the radial line profiles if the azimuthal line profiles are 
given. This procedure is illustrated on [Fig. 2.3 for simulated data without noise. The surface density £ 
follows from both Fr(vr, R) and F^v^, R) and could also be compared with the photometry of that galaxy. 
The likely discrepancy between the predicted and the residual data may be used to asses the limitations of 
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the reconstruction scheme. Indeed, the above prescription for the inversion relies on a set of hypotheses for 
the nature of the flow (axial symmetry, thin disk approximation etc..) and the quality of measurements. 
The following features of an observed galaxy may constrain the scope of this analysis. 
• Thick disk, with random motions perpendicular to the plane of the galaxy. 

The measured line profile yields a mean value corresponding to the integrated emission through the 
width of the disk; as the galaxy is not edge on, the result of the cross correlation of the spectra does not 
give Fcj, and Fr exactly, but rather gives the projection onto the line of sight of the full velocity distributions 
-p0,z(w0, v z , R) and Fr iZ (vr, v z ,R). However, it is consistent with the thin disk approximation to assume that 
the motion in the plane is decoupled from that perpendicular to the plane. The observed line profile, F v [R, v), 
then corresponds to the convolution of these functions with the component of the velocity distribution 
perpendicular to the plane of the galaxy, F Z (R, v z ). Formally, for the line profile measured along the major 
axis, this reads 



where v and u are the velocities along the line of sight, and perpendicular to the line of sight respectively. 
A similar expression for the line profile measured along the minor axis involving Fr follows. Here, the 
angle i measures the inclination of the plane of the galaxy with respect to the line of sight. At this level of 
approximation the velocity distribution normal to the plane is accurately described by a centered Gaussian 
distribution; its variance follows roughly from the equation of radial support and is used to deconvolve 
Eq. (4.1) (or practically to convolve the predicted line profile pairs and compare them to the data). Note 
that a rough estimate of the error induced is ct^/ct? tan 2 (i) ~ 1% in our Galaxy when i = 20 degrees. Note 
that when looking at an ellipse that corresponds to a circle in the galaxy's plane one would measure the 
dispersion (cfj + cr 2 ,) cos 2 (z)/2 + (cr R — cr?) cos(2</>) cos 2 (i)/2 + <r 2 sin 2 (z). If it is assumed that the asymmetric 
drift hypothesis holds, - which can be checked by measuring (v^) - then it should in principle be possible 
to measure indirectly a z . 

• Non-axisymmetric halo or disk, with dust and absorption. 

The departure from axisymmetry reduces the number of invariants - angular momentum is not con- 
served, which invalidates this method. Some fraction of the population of stars may then follow chaotic 
orbits which cannot be characterised by a distribution function. Some real galaxies either present intrinsic 
non axisymmetric features such as bars, spiral or lopsided arms. Triaxial halos may induce warps within 
the disk, hence compromising the evaluation of the kinematics. Non-axisymmetrically distributed dust or 
molecular clouds may also affect measurement of the light from the stellar disk. On the positive side, non- 
axisymmetric features often appear as small perturbations which may only contribute weakly to the overall 
mass distribution. Under such circumstances, an estimate of the underlying distribution function for the 
axisymmetric component of the disk may be extracted from the data. 

Sets of initial conditions for an N-body simulation may also be extracted from the distribution function. 
The time evolution of the code would yield constraints on the stability of the observed galaxy. The above 
analysis could be carried out on an isolated SO galaxy and on an Sb galaxy showing a grand design spiral 
pattern (or alternatively on an Sa presenting an apparent lopsided HI component to isolate plausible different 
instability modes). The relative properties of the two galaxies and their supposedly distinct fates when 
evolved forward in time should provide an unprecedented link between the theory of galactic disks and the 
detailed observation of the kinematics of these objects. The authors fl5| have written a numerical procedure 
to implement linear stability analysis on arbitrary mass models and kinematics for flat and round galactic 
disks. 

• Data quality. 

The method described in section 2 involves a deconvolution of the measured line profiles which yields the 
velocity distributions which, in turn, are related by Abel transforms. Both transformations are very sensitive 
to the unavoidable noise in the data. Sources of noise are poor seeing, photon noise, detector noise, thermal 
and mechanical drift of the spectrograph, poor telescope tracking, etc... It is therefore desirable to construct 
from Eqs. (2.4), (2.6) and (2.17) a set of pairs of Gaussian velocity distributions (F^, Fr) on which to project 
the observed quantities since this will help to assess the errors induced by this noise. Indeed, Eq. (2.4) is 
linear; therefore any superposition of Gaussian profiles corresponding to a good fit of the observed line profile 




F z [R, v sin (i) + u cos (i)] F^ [R, v cos (i) — u sin (i)] du 



(4.1) 
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will lead to a distribution function expressed in terms of sums of solutions of Eq. (2.4). These Gaussians may 
in turn be identified with populations at different temperatures - say corresponding to a young component 
and an older population of stars. Formally this translates as 



where the functions Si(R),(ii(R) and Vi(R) should be fitted to the observed azimuthal line profiles. Each 
component in Eq. (4.2) may then be identified with those in Eq. (2.20), yielding the temperature profile 
of that population. This approach should yield the best compromise between fitting the radial and the 
azimuthal profile while avoiding the deconvolution of noisy data. It also provides directly a least squares 
estimate of the errors. The observations described in this section may be achieved with today's technology. 

5. Conclusions 

Simple general inversion formulae to construct distribution functions for flat systems whose surface density 
and Toomre's Q number profiles were given and illustrated. The purpose of these functions is to provide 
plausible galactic models and assess their critical stability with respect to global non-axisymmetric pertur- 
bations. The inversion is carried out for a given azimuthal velocity distribution (or a given specific energy 
distribution) which may either be observed or chosen accordingly. When the azimuthal velocity distribution 
is measured from data, only a subset of the observationally available line profiles, namely the line profiles 
measured on the major axis, is required to re-derive the complete kinematics from which the line profile 
measured along the minor axis is predicted. This prediction may then be compared with the observed minor 
axis line profile. The Ansatz presented in section 3 yield direct and general methods for the construction 
of distribution functions compatible with a given surface density for the purpose of theoretical modelling. 
These distribution functions describe stable models with realistic velocity distributions for power law disks, 
and also the Isochrone and the Kuzmin disks. Note that the inversion technique described in section 2 can 
be implemented in the relativistic regime as discussed by Pichon & Lynden-Bell ]l6| . 

CP wishes to thank R. Cannon and A. Cooke for usefull conversations, and the anonymous referee for 
suggesting improvements in the organisation of the paper 
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